Evaluation of Exchange-Correlation Energy, Potential, and Stress 
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We describe a method for calculating the exchange and correlation (XC) contributions to the 
total energy, effective potential, and stress tensor in the generalized gradient approximation. We 
avoid using the analytical expressions for the functional derivatives of _B xc (p), which depend on 
discontinuous second-order derivatives of the electron density p. Instead, we first approximate _B XC 
by its integral in a real space grid, and then we evaluate its partial derivatives with respect to 
the density at the grid points. This ensures the exact consistency between the calculated total 
energy, potential, and stress, and it avoids the need of second-order derivatives. We show a few 
applications of the method, which requires only the value of the (spin) electron density in a grid 
(possibly nonuniform) and returns a conventional (local) XC potential. 
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I. INTRODUCTION 

The generalized gradient approximation (GGA)§ to 
density functional theory (DFT)B has been growing in 
acceptance in recent years, due to the development of 
improved functionals and to the realization of its higher 
accuracy, for many systems and properties, than that of 
the local density approximation (LEjA)EI. Some continu- 
ity problems of earlier functionalstffl have been solved 
in more recent onesu. Still, a basic problem is that, 
while E t £p A [p] depends only on the local value of p(r) 
and |Vp(r)|, its functional derivative depends also on 
V|Vp(r)|, which is discontinuous where Vp(r) = 0. This 
implies that, to avoid aliasing effects, very fine integra- 
tion grids are required to evaluate the XC potential u xc (r) 
and its matrix elements. A different problem occurs with 
grid-oriented implementations of DFTmij, in which the 
electron density is known only at the grid points, while 
its gradients must be evaluated using finite differences. 
In this case, an inconsistency between the energy E xc 
and the potential v xc may occur because different formu- 
las need to be used for the higher-order derivatives of the 
density. Hamannu proposed an elegant solution to these 
problems by defining a nonlocal potential, which operates 
on the gradient of the electron wave functions. Although 
in principle this does not pose any special difficulty, it re- 
quires a specific, unconventional implementation, which 
may be difficult to adapt to existing DFT codes. White 
and Birdt2l found another solution, by defining E^F A as 
an integral in a real space grid (what is always done in 
practice), and v xc (i"i) as its partial derivative (as opposed 
to functional derivative) with respect to the density pi at 
the grid points rj. They applied this method to a plane 
wave basis set and a uniform integration grid, finding the 



gradients with fast Fourier transforms (FFTs). Here we 
generalize their method to arbitrary bases and nonuni- 
form grids, and we calculate the density gradient using 
finite differences. The method produces a standard, lo- 
cal potential, which is exactly consistent with the defini- 
tion of -E xc , in the sense that v KC is the correct potential 
in the Schrodinger equation that results from the varia- 
tional minimization of the total energy. In addition, we 
show how to evaluate the XC contribution to the stress 
tensor in crystals, following the same ideas. The trivial 
particular case of the LDA is also discussed. 



II. EXCHANGE-CORRELATION ENERGY AND 
POTENTIAL 

The LDa! and GGA0 approximations to E xc have the 
general forms 



E^ DA [P}= JhDA(p(r))dr, 
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/GGA(p(r),g(r))dr, 



(1) 



where / = p e xc is the local XC energy density, and f__ 
is the XC energy per electron. We use the notationli^l 
g(r) = Vp(r), g a (r) = V Q p(r), and g(r) = |Vp(r)|. The 
argument g(r) in Eq.(|l]), rather than g(r), indicates that, 
in principle, /gga might depend on the relative orien- 
tation of the gradients of the two spin components of 
the density. However, for notational simplicity we will 
omit the spin index and the sum over it, since their 
inclusion is trivial. In the limit of slowly varying elec- 
tron density, we must recover the LDA result, that is 
/gga (p, 0) = /ldaG°)H. 
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In the LDA the XC potential has a simple expression, 
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whereas in the GGA it is considerably more complicated 
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as it depends on the first and second gradients of the elec- 
tron density. Since V/?| has cusps at extrema of p, Vg is 
discontinuous at these points, what causes problems for 
its numerical representation. Furthermore, some param- 
eterizations of 6xc GA which apparently join seamlessly dif- 
ferent density regimes, may have higher derivatives which 
do not behave well in those transition regions. 

In what follows, we will omit the label GGA, except to 
underline a distinction with LDA. In a practical calcu- 
lation, the XC energy is calculated through a numerical 
integration. From a set of M mesh points and weights 
Wi, we approximate 



M 

£ 

i=l 



Wi/(p(r. ( ),g(r 4 )). 



The weights may be, for example AirrfAr in a uniform 
radial grid, or the jacobian .a£ the local metric tensor in 
an adaptive-coordinate gricO (see below). In addition, 
we must specify precisely the meaning of g(i"j). If a well 
defined basis set is used, the electron density and its gra- 
dient can be calculated exactly at any point in space, 
from the electron wave functions and their gradients. In 
practice, this may add an appreciable overhead in terms 
of computer time and memory. Alternatively, we can use 
the values of the density at the, grid points to calculate 
its gradient, using either FFTsEJ or finite differences. As 
the gradient is a linear operator, we can write in general 
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and dgi a /dpj — Afj. {i,j} are combined indexes (i = 
{ill *2> *3}) that label grid points, a labels the three carte- 
sian coordinates, pi = p^i), and the coefficients Afj de- 
pend on the mesh and the chosen numerical derivative 
formula, but not on pi. There are many choices for the 
coefficients Afp depending on the integration mesh and 
interpolation method, since the numerical derivative can 
be defined as the derivative of the interpolation function. 



The particular case described by White and BirdEl uses a 
uniform mesh and a Fourier-series interpolation. It was 
developed for plane-wave calculations, and FFTs were 
used to evaluate Eq. (||) . Our method can be generalized 
to nonuniform grids. In the implementation that we will 
describe later, we use a local Lagrange interpolation, for 
which the matrix Afj is sparse. 

Notice that in Eq. (§) we are defining gi as the nu- 
merical derivative of p on the mesh (therefore being a 
function of the values {pj}), while we reserve the nota- 
tion g(rj) for the exact gradient of p at (in case it is 
known). It may be argued that the use of gi, instead 
of g(i"j) represents an additional approximation. How- 
ever, in practice gi is frequently a better approximation 
than g(rj) to the average value of g(r) within the spatial 
"pixel" which corresponds to the integration point rj. 
Both values agree for a Fourier interpolation, provided 
that the plane wave cutoff of the grid is twice as large 
as that of the wave functions. And of course they must 
also agree, for any interpolation scheme, in the limit of 
an infinitely fine grid. _ . 

Thus, following White and Birdli 2 ], we define 
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and 
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It is important to emphasize that Eqs. ^ and (||) are 
asymptotically equivalent in the limit of an infinitely fine 
grid, but different in practice (in particular, Eq. (^) has 
no explicit dependence on the second derivatives of p). 
In fact, it is easy to see that Eq. (||), and not Eq. (j^), is 
the "correct" definition of v xc , if the functional form (Q) 
is actually used in the variational energy minimization, 
dE/difj* = Hi/j, because, from Eq. (|j), 



dE xc _ dE xc dp l 
dip* dpi dipi 



In the LDA, /; does not depend on gi, and only the 
first term in Eq. (||) remains, giving the trivial result 
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There remains the problem of determining the coeffi- 
cients Afj and the weights Wi in Eqs. (||) and (^). To 
this purpose, it is convenient to introduce grid vari- 
ables {s^p — 1,2,3}, which are in principle contin- 
uous. In practice, however, the density p(s) and the 
cartesian coordinates r(s) of the grid points are evalu- 
ated only at integer values of s M . For a regular grid, 
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r a( s ) — J2 t i=i s ^ a tJ-a/N fl , where a MCt is the ath carte- 
sian coordinate of the pth lattice unit vector, and 
is the number of grid divisions along that vector. For 
nonuniform grids, r Q (s) are general functions and it is 
convenient to introduce the jacobian of the transforma- 
tion s — > r, The XC energy can then be expressed as 



E, c = J /(p(s),g(s)) 
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or 
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point j: 1) find the 3x3 matrix Dj, using Eq. (^|), its 
determinant Wj , and its inverse DJ , storing Wj for later 
use; 2) find A£ from Eq. @) and gj from Eq. (g); 3) 
calculate f(pj,gj) an d its derivatives with respect to /X; 
and g ja ; 4) add the first term in Eq. d|) (multiplied by 
Wj) to and the second term (except the denominator 
w%) to -Df c for every neighbor point i involved in the cal- 
culation of gj; 5) when the previous loop is finished, run 
again over all grid points i, dividing vf c by Wi. In the 
case of a uniform grid, the matrix Di does not depend 
on i. Thus, wt and D~ x can be evaluated once and for 
all, saving steps 1 and 5, as well as the temporary array 
required to store w;. 



where we have used that As M = 1 by definition of s M . 
|£>i| is the determinant at point i of the matrix of partial 
derivatives 



D aiJ ,(s) 



dr a (s) 
ds,. 



By comparison with (0), we conclude that Wi = |-D,|. 
Since the grid variables s^ are by definition "orthogonal" 
and the functions r a and p are evaluated at regular unit 
intervals of s M , their derivatives d/ds^ can be calculated 
straightforwardly: 



III. EXCHANGE-CORRELATION 
CONTRIBUTION TO THE STRESS TENSOR 

We consider now the stress tensoill! 

dE xc 



'a/3 



de a 



13 



where e a p is the strain tensor, giving the deformation 
of all points in space (including atomic and grid coordi- 
nates): 
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ds u 



M 

4=1 



(6) 



0=i 



(8) 



ds. 



A I 



3 = 1 



where the coefficients B^ depend only on the relative 
values of i and j, and are independent of the grid coor- 
dinates r(s): 
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where the coefficients may be derived from a (JZn+l) 
point Lagrange polynomial interpolation formulaEa. 
We can now calculate g as 
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Thus, the calculation of the XC potential from the den- 
sity on the grid involves the following steps, for every 



More generally, we consider the derivative of E xc with 
respect to a parameter A that affects the position of the 
grid points. This may be e a p or one of the atomic po- 
sition*, in case of a nonuniform grid which depends on 
themlU. It is therefore convenient to recognize explicitly 
that, when the system is modified, E xc depends on the 
grid point coordinates in addition to the densities at the 
grid points, i.e. E xc ({pi}, {i-j}). Thus, 
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Now, using Eq. QSj) we find 
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and taking into account the general properties of a deter- 
minant, as well as the fact that (DD^ 1 = I) =>■ (5D = 
-D~ 1 8DD~ 1 ), we finally find 
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In the limit of an infinitely fine grid, we have 
= S a pE KC + / v xc {r)— dr 



de af: 
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The presence of dp/e a p in Eqs. (Jlfj) and (|llj) requires 
to clarify what is kept constant when taking the deriva- 
tives. We will argue that the "correct" definition (or 
at least the most useful one) requires to keep constant 
the variational coefficients c na of the expansion of the 
electron wave functions ip n (r) in terms of the N basis 
functions 4>a( r )' 
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ip n {r) = ^c„ a Q (r) 



(12) 



Since a a p is generally evaluated at the electronic ground 
state, the Hellman-Feynman theorem ensures that the 
change of c na will not affect the total energy to first order. 
Of course, the Hellman-Feynman theorem does not apply 
to E xc alone, but in practice we are interested in adding 
a af3 t° other contributions of cr a p, in order to calculate 
the total stress tensor. And even in the Car-Parrinelldl3 
method, in which forces and stresses are evaluated out 
of the ground state, the appropriate definition of these 
magnitudes involves the derivatives of the total energy 
at constant c na . 

There are two reasons why pi depends on t a p: the 
change of the basis functions </> (rj) at the displaced grid 
points, and the change of the wave function coefficients 
Cna required to maintain the orthonormality constraints 
{ipn\4>m) — finm- With,.-a-jplane wave basis set, or with 
a grid-oriented schemeoEljiJ, in which the wave functions 
arc defined directly at the grid points, without any ba- 
sis, the orthonormality constraints are not affected by 
the deformation of the unit cell, provided that the wave 
functions and the density are simply scaled by a factor 
(Qq/Q) 1 / 2 and £Iq/SI, respectively, where f2o and Q are 
the unit cell volumes before and after the deformation. 
It is easy to verify that dQ/de a p = Vl8 a p and therefore 
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Substitution of ( |l3| ) into (^lj) leads to 
dE xc 
de a /3 



= S af 3 / (e xc (r) - -u xc (r))p(rrfr 
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which coincides with Eq. (24) of Dal Corso and RestaEl. 

In the case of an atomic basis set, it is convenient to 
define the density matrix 



N 



Pab ^ Qn^na^nb 
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with q n the occupation of state ip n . Then 
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P(r) = E /W>a(r)0 fc (r) 
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where we are assuming real basis orbitals for simplicity. 
Rb is the origin (atomic position) of orbital (f> b , and the 
last term in Eq. (jlj) accounts for the change in the rela- 
tive position r — Rj, when we deform the lattice or change 
the grid point positions (we assume a constant shape of 
4>b)- dpab/dX is the change in the density matrix required 
to maintain the orthonormality constraints: 
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where S ab = (4> a \4>b), and r ab = R b - R a . By combining 
Eqs. @, @, and @ we obtain 

v xc {r)- dr = - > p ac -K-^r>; b (4> a \v xc \<j> b ) 

J oe a p f-^ dr% 

a,6,c=l ' 
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The same expression is valid if the integral is replaced by 
the grid sum J2i w ivf°{dpi/ d^af}), except that then the 
matrix elements must also be calculated on the grid, i.e. 
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f>a\v xc \(f>b) = ^ WjVf C (j) al 



i=l 

Similar terms appear in other contributions to the stress 
tensor. Thus, in the derivation of the Hartree energy, 
there is an identical term, except for the substitution of 
^xc by the Hartree potential vh- Then the total contri- 
bution to the stress tensor that arises from the change in 
the density has the form 



Sp(r) de af 3 
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where v(r) is the total effective potential, and E a b are 
energy-density matrix elements, 
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E a b = ^ Pa 

c=l 



H c b = 
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(17) 



with H a b — (<t>a\H\ < t ) b}, H the total one-electron hamil- 
tonian and e„ its nth eigenvalue. Notice, however, that 
the derivation of Eqs. ( [l6| ) and ( O ) does not assume that 
the matrix elements S^f, and H a b or their derivatives are 
calculated on tlie,grid. In fact, in our atomic-basis DFT 
implementationE2l, two-center matrix elements like S a b or 
((j>a \ — V 2 \(j)b), are calculated using reciprocal-space con- 
volution techniques^. 



IV. ACCURACY TESTS 

We have implemented the method described above as 
an independent package, which we are using with sev- 
eral electronic structure programs, including an atomic 
pseudopptential generator, the localized-basis program 
SiestaEJ and a pseudopotential plane-wave program. 

We have implemented two different LDA parameter- 
izations, PZ-LDAE3 and PW92-LDAy, of Ceperley and 



Adler electron gas energies 
GGA recipes, PW91-GG 



cP3| £ or j LDA an d Hir different 
and PBE-GGAQ for / GGA . 
They were all implemented in their spin-dependent for- 
mulations. As the calculation of / and its derivatives are 
kept as separate procedures it is easy to add new param- 
eterizations. As only first derivatives of / are needed in 
our method to calculate stresses and potentials, we do 
not have to calculate second derivatives as in the tradi- 
tional method. Since GGA functional forms tend to be 
complicated, this is convenient. 

In the atomic pseudopotential generation program the 
electron density is assumed to have spherical symmetry. 
The radial mesh points are derived from a monotonous 
function r(s). The gradient of the density is calculated 
analogously to the three-dimensional case: 
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To illustrate the method we have chosen three sim- 
ple situations. Fig. [I] shows the s-jaseudopotential for 
Si, with the Troullier-Martins recipea and a core radius 
of two bohr, calculated with the PBE-GGA and PW91- 
GGA functionals using both our method (Eq. (j^)) and 




r (a.u.) 

FIG. 1. s-pseudopotential generated for Si. The solid 
curve was generated with the PW91-GGA and the dashed 
curve with the PBE-GGA. The curves generated with the 
traditional and the new methods are undistinguishable. 



the usual method (Eq. (g)), with the standard radial grid 
used by the pseudopotential generation program. The 
curves for the two computational methods agree within 
10~ 4 Ry and are indistinguishable in the scale of the fig- 
ure. The PW91-GGA pseudopotential has some wiggles 
that are due to the term proportional to exp(— 100s 2 ) 
where s = Vp/(2fcpp), fcp = (S^p) 1 ^ 3 , in the parame- 
terization of the GGA exchange. Those wiggles occur 
near the extrema of the valence electron density, and 
they are not present in the PBE-GGA pseudopotential, 
which uses a different parameterization of the GGA ex- 
change. By construction the screened pseudopotential 
is smooth, without wiggles, so the appearance of those 
wiggles are due to the unscreening of the pseudopoten- 
tial. To see any difference between the usual and present 
methods, we must use very coarse radial meshes. This is 
shown in Fig. |^, where we compare the XC potential for 
a model density p(r) — (sinr/r) 2 . The PZ-LDA (dotted 
thin line), PW91-GGA (solid thick line), and PBE-GGA 
(dashed line), were first calculated on a fine radial mesh. 
Again there is no difference between the old and new 
methods on the scale of fig. || for the fine mesh. How- 
ever, recalculating the PW91-GGA exchange and corre- 
lation potential on a much coarser mesh with the new 
method (squares) one starts to see small differences in 
the potential values in the regions where the PW91 pa- 
rameterization has wiggles. The eleven point Lagrange 
formula for coarse sampling spans to a region of 2 bohr, 
almost half of the horizontal range of Fig. 0. It misses 
the small wiggles around r — 3tt/2 where Vp = and 
p 7^ 0, but it finds wiggles at r = ir where both Vp = 
and p = 0. 

Our last example uses the pseudo-charge density of di- 
amond. A well converged plane wave expansion of the 
density is used to calculate pi for different N x N x N 
grids. Fig. || shows that the calculated XC energy is 
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FIG. 2. XC potential for a spherical electron density with 
the form p(r) = (sin r/r) . The dashed line is the CA-LDA, 
the dotted line is the PBE-GGA and the solid line is the 
PW91-GGA, all accurately calculated on a fine grid with a 
step of 0.02 bohr. The squares show that the PW91-GGA, 
calculated with the new method on a coarse grid (step = 
0.2 bohr) have small deviations with respect to the fine grid. 
However, in fact these deviations tend to smooth out the large 
pathological wiggles developed by the PW91-GGA functional 
in regions of zero density or zero density gradient. This makes 
the present method very well behaved even with very coarse 
grids. 
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FIG. 3. XC energy in diamond as a function of the grid 
size N x N x N. Squares are from LDA and circles from 
GGA. The charge density p(r) is the same in all cases. Only 
the sampling grid and the XC functional change. 
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very stable with respect to grid size: even a 8 x 8 x 8 
grid gives a value accurate within 10 -3 Hartree. The XC 
component of the stress also converges very fast, and the 
method is very stable. In practice, the PW cutoffs re- 
quired for a good convergence of the total energy impose 
larger grids than those needed to converge E xc with our 
method. Thus, the standard grids used in typical PW 
calculations are more than sufficient for E xc . 



V. CONCLUSIONS 

We have derived formulas for the GGA-XC energy, po- 
tential, and stress, as a function of the electron density 
given in a spatial mesh of points, which may be unevenly 
distributed. Density gradients need not be provided, 
since they are calculated numerically using the density 
at the grid points, what leads to well behaved formulas. 
The use of a unique definition of the gradients ensures an 
exact consistency between the calculated values for the 
energy, potential, and stress, including all corrections for 
basis set incompleteness and changing grid points. As 
the number of grid points increases, we recover the exact 
results obtained using the virial theorem. 



1 D. C. Langreth and M. J. Mehl, Phys. Rev. B 28, 1809 
(1983). A. D. Becke, Phys. Rev. A, 38, 3098 (1988). 

2 P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964). 

3 W. Kohn and L. J. Sham, Phys. Rev. 136, A1133 (1965). 

4 J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). 
J. P. Perdew in Electronic Structure of Solids 91, p. 11, 
edited by P. Ziesche and H. Eschrig (Akademie Verlag, 
Berlin, 1991). 

5 C. Filippi, C. J. Umrigar and M. Taut, J. Chem. Phys. 
100, 1290 (1994). C. Filippi, X. Gonzc and C. J. Umrigar 
in Recent Developments and Applications of Modern Den- 
sity Functional Theory, p. 295, edited by J. M. Seminario 
(Elsevier, 1996). 

6 J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 
77, 3865 (1996). We are not considering a new generation of 
GGA funtionals including a dependence with the laplacian 
of the electron density. Nevertheless, the main approach in 
this paper is even more useful for that type of functionals. 
See, i.e., A. Rassolov, J. A. Pople and M. A. Ratner, Phys. 
Rev. B 59, 15625 (1999); ibid 62, 2232 (2000). 

7 J. R. Chelikowsky, N. Troullier, and Y. Saad, Phys. Rev. 
Lett. 72, 1240 (1994); J. R. Chelikowsky, N. Troullier, K. 
Wu, and Y. Saad, Phys. Rev. B 50, 11355 (1994). 

8 D. R. Hamann, Phys. Rev. B 54, 1568 (1996). 

9 N. A. Modine, G. Zumbach, and E. Kaxiras, Phys. Rev. B 
55, 10289 (1997). 



G 



10 E. L. Briggs, D. J. Sullivan, and J. Bernholc, Phys. Rev. B 
52, R5471 (1995); ibid 54, 14362 (1996). 

11 T. L. Beck, Rev. Modern Phys. 72, 1041 (2000) 

12 MJ. A. White and D. M. Bird, Phys. Rev. B 50, 4954 
(1994). 

13 For convenience, we will place indexes and labels indis- 
tinctly as sub or superscripts, without meaning a tensorial 
notation nor implicit summations. Indexes i, j are used for 
the M mesh points, a, f3, 7 for the three cartesian coordi- 
nates, fi, v for the three lattice or mesh coordinates, a, b, c 
for the N basis functions, and n, m for the N hamiltonian 
eigenstates. 

14 F. Gygi,Phys. Rev. B 48, 11692 (1993); F. Gygi and G. 
Galli, Phys. Rev. B 52, R2229 (1995). 

15 We use here the convention that the stress has the opposite 
sign than the pressure. Often the opposite convention is 
used, so the reader may find an overall change of sign in 
the expressions given by different authors. 

16 As the energy of a crystal does not depend on its orienta- 
tion in space, only the length of the lattice vectors and the 
angles between them are important. In other words, the 
metric tensor G M „ = a M • a„ contains all the relevant in- 
formation. This allows an equivalent formulation, in which 
the contravariant components of the stress tensor can be 
denned as r MI/ = 2dE / 'dG M „ and its cartesian components 
are then a a p = ^^T^a^aCivp. For details, see I. Souza 
and J. L. Martins, Phys. Rev. B 55, 8733 (1997). 

17 O. H. Nielsen and R. M. Martin, Phys. Rev. B 32, 3780 
(1985). 

18 A. Dal Corso and R. Resta, Phys. Rev. B 50, 4327 (1994). 

19 R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985). 

20 P. Ordejon, E. Artacho, and J. M. Soler, Phys. Rev. B 53, 
R10441 (1996); D. Sanchez-Portal, P. Ordejon, E. Artacho, 
and J. M. Soler, Int. J. Quant. Chem. 65, 453 (1997). 

21 O. F. Sankey and D. J. Niklcwski , Phys. Rev. B 40, 3979 
(1989). 

22 J. P. Perdew and A. Zunger, Phys. Rev. B 23 , 5048 (1981). 

23 D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45 , 566 
(1980). 

24 Handbook of Mathematical Tables, Ed. by M. Abramowitz 
and I. A. Stegun (Dover, New York, 1965). 

25 N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 
(1991). 



7 



